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^T*' By means of molecular dynamics computer simulations we investigate the out of equilibrium 

relaxation dynamics of a simple glass former, a binary Lennard- Jones system, after a quench to 
low temperatures. We find that one-time quantities, such as the energy or the structure factor, 
show only a weaJi time dependence. By comparing the out of equilibrium structure factor with 
equilibrium data we find evidence that during the aging process the system remains in that part of 
phase space that mode-coupling theory classifies as liquid like. Two-times correlation functions show 
a strong time and waiting time dependence. For large and times corresponding to the early 
^ ' /3-relaxation regime the correlators approach the Edwards- Anderson value by means of a power-law 

Ch . in time. For large but fixed values of the relaxation dynamics in the /3-relaxation regime seems to 

^ ' be independent of the observable and temperature. The a-relaxation shows a power-law dependence 

on time with an exponent which is independent of but depends on the observable. We find that 
at long times r the correlation functions can be expressed as C'AG{h{tw + t) / h{tw)) and compute the 
function h{t). This function is found to show a t-dependence which is a bit stronger than a logarithm 
and to depend on the observable considered. If the system is quenched to very low temperatures the 
relaxation dynamics at long times shows fast drops as a function of time. We relate these drops to 
I relatively local rearangements in which part of the sample relaxes its stress by a collective motion 

(~j , of 50-100 particles. Finally we discuss our measurements of the time dependent response function. 

Q ' We find that at long times the correlation functions and the response are not related by the usual 

O , fluctuation dissipation theorem but that this relation is similar to the one found for spin glasses 

with one step replica symmetry breaking. 
PACS numbers: 61.20.Lc, 61.20.Ja, 02.70.Ns, 64.70.Pf 
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'sj- ' I. INTRODUCTION 

(N ■ 

, It is well known from experiments that the most prominent feature of supercooled liquids and glass forming systems 
' is the rapid increase of their relaxation time as temperature is decreased Since the early eighties, computer 
^\ , simulations have become a powerful tool to investigate such systems and hence to increase our understanding of their 
extraordinary dynamical behavior Q. However, because of the rapid increase of the relaxation times these simulations 
are faced with the intrinsic problem that, over a large range of temperature, systems cannot be equilibrated, as the 
simulation time, typically 10~^s for an atomic system, is smaller than the relaxation time. This means that on 
I a computer, any system with a longer relaxation time can only be studied under nonequilibrium conditions. In 
other words, in computer simulations the "glass transition" always takes place at temperatures higher than in real 
experiments. 

^ , Despite this drawback, it has been shown that computer simulations can provide useful information on the mccha- 
' nisms than underly the rapid increase of the time scales in supercooled liquids. One of the most valuable contributions 
. !^ ' in this direction was the proof that the so called mode-coupling theory does indeed provide a quantitative de- 
scription of this slowing down within the time window explored in the simulation, at least for those simple atomic 
liquids that are "fragile glass formers" , i.e. show a strongly non-Arrhenius dependence of the relaxation times as a 
function of temperature. The findings from computer simulations ||^ and experiments Q have largely confirmed that 
this theory provides a consistent picture for the dynamical behaviour of such liquids over a relaxation time range that 
covers several decades, typically 10~^^ — 10~®s. 

For larger relaxation times (lower temperatures), equilibrium simulations are no longer possible. On the other 
hand, one can take advantage of the flexibility of the simulation to explore the nonequilibrium properties on times 
scales that would be rather difficult to access in experiments. This possibility stems from the fact that the simulations 
allows instantaneous quenches to low temperatures, and a monitoring of the subsequent evolution on a fast time scale. 
Experimentally, the study of nonequilibrium evolution in glassy systems following a fast quench is a field that has been 
explored since a long time in the so called "aging" experiments pi, performed on a time scale of hours or even days. 
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More specifically, "physical aging" consists in a slow evolution of the characteristic properties of glassy systems that 
do not undergo any changes in their chemical composition. Of particular interest is the evolution of response functions 
(typically the elastic compliance) that are usually found, e.g. in polymeric systems, to exhibit a strong slowing down 
as the waiting time (i.e. the time elapsed since the quench) increases. The recent interest of the physics community 
in this behavior was aroused when a number of experiments performed on spin glasses showed that a very similar 
behavior could be observed in these disordered magnetic systems. In particular it was found that the response to a 
magnetic field of a system that has been quenched into its low temperature, nonequilibrium spin glass phase, becomes 
more and more sluggish as the waiting time increases Subsequently, a number of theories were put forward to 
explain this aging behavior of spin glass systems from first principles statistical physics. A comprehensive review of 
these theories can be found in reference 0. Essentially, one is lead to distinguish between phenomenological "trap 
models", domain growth theories, and mode-coupling or mean field theories. "Trap models" describe the evolution of 
systems as a random walk in a complex phase space, which can under certain conditions give rise to aging. The domain 
growth models assume that the aging results from a coarsening process somewhat similar to what can be observed in a 
ferromagnetic system quenched below its critical temperature. Mean field or mode-coupling theories account exactly 
for the aging behaviour of some disordered models in the limit of high dimension (e.g. a particle in a random potential 
in a high dimensional space). Interestingly, these mean field theories give rise to a mathematical structure that is 
very similar to that of the mode-coupling theory commonly used to describe the dynamical properties of supercooled 
liquids A major difference, however, is that high dimensional systems can be shown to have a true transition 

towards a non-ergodic behavior as the temperature is lowered, whereas it is well known that the transition predicted 
by the mode-coupling theory of liquids is smeared out and does not give rise to a true singularity. If, however, a 
true transition is present, a natural consequence is that a system quenched below the transition temperature will 
display aging behavior. Theoretically, this aging behavior of mean field systems was studied in great detail and is 
now quite well understood [Q. The question arises, however, to what extent this description is relevant to the aging 
behaviour of real, three dimensional systems. In particular, in view of the mathematical similarity mentioned above, 
it is quite natural to inquire whether the nonequilibrium dynamics of those supercooled liquids that are known to be 
well described by mode-coupling theory can be accounted for by the mean field description of aging. 

In an attempt to partially answer this question, we present molecular dynamics simulations of a simple glass 
forming liquid, whose properties under equilibrium conditions have been extensively studied in previous work 
We will concentrate on the characterization of the "aging" behaviour of such a system on the time scales that can 
be investigated in molecular dynamics simulations, i.e. a few nanoseconds. Once again, we emphasize that these 
time scales are very different from those usually investigated in aging experiments on glasses (although some recent 
dielectric spectroscopy experiments [p|JTo|] investigate relatively high frequency behaviour). Moreover, the systems 
under study would, from an experimental viewpoint, be considered as supercooled liquids rather than glasses. Glassy 
behavior is observed only on the restricted time scale of the computer simulation. The possibility of exploring such 
relatively short time scales is nevertheless interesting, as it is precisely on such time scales that mode-coupling theory 
is successful. Moreover, the simulation offers the opportunity to simultaneously compute time dependent correlation 
functions, static quantities and response functions. Relating these quantities to each other is an essential achievement 
in mean field theories of aging behaviour, whose predictions can then be tested in detail. 

The paper is organized as follows. In section II, we recall the main features of our model and describe some 
technical details of the simulations. Sections III and IV describe the aging behavior of static properties and of 
correlation functions, respectively. Section V deals with the aging behaviour of response functions. Preliminary 
accounts of these results can be found in Some closely related studies, on a slightly different system, have been 
reported in 



II. MODEL AND DETAILS OF THE SIMULATION 



The model we use is a binary (80:20) mixture of particles, which in the following we will call A and B particles, 
that interact via a Lennard- Jones potential, Vaf){r) — ^tapii'^ap/TY'^ — [^ap/f')^}^ with a, /3 g {A, B}. The constants 
tap and UajS are given by eAA = 1-0, ctaa ~ 1-0, cab — 1-5, ctab — 0.8, eeB — 0.5, and ctbb = 0.88 and the potential 
is cut off and shifted at a distance 2.5tTQ/3. In the following we will report the results in reduced units, with ctaa and 
EAA the unit of length and energy, respectively (setting the Boltzmann constant ks — 1.0). Time will be measured in 
units of cr^^m/48eAA, where m is the mass of the particles. 

The simulation was done with N = 1000 particles and at constant volume V = L^, with a box length L = 9.4. 
At this density the equilibrium dynamics of the system has been investigated intensively [^|JT^ and it has been found 
that at low temperatures, 0.446 < T < 0.8, this dynamics is described very well by mode-coupling theory [^||^] 
with a critical temperature Tc ~ 0.435. Although there is evidence that for temperatures T < 0.452 the increase of 
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the relaxation times, or the inverse of the diffusion constant, is less strong than the power-law predicted by MCT, 
(T — Tc)~'^ , with 7 = 2.3 the increase is still so strong that at temperatures less than Tc the relaxation time of 
the system exceeds by orders of magnitude the time scale accessible to a present state of the art computer simulation 
(0(10'') time units). Hence for a computer simulation Tc plays the role of the glass transition temperature Tg. 

The equations of motion were integrated with the velocity form of the Verlet algorithm, with a step size of 0.02. 
The system was equilibrated in the high temperature phase at a temperature Ti > Tc and then quenched at time 
t = to a temperature T/ < Tc- This quench was done by means of a stochastic heat bath, i.e. every 50 time 
steps we replaced the velocities of all the particles with ones that were drawn from a Maxwell-Boltzmann distribution 
corresponding to a temperature T/ . In order to study the dependence of the relaxation behavior on the initial and 
final temperature we considered several values of Ti and Tf. Ti £ {5.0,0.8,0.466} and Tf £ {0.435,0.4,0.3,0.1}. 

After the quench we propagated the system for a waiting time ty^ after which the measurements of the quantities 
of interest were started. In the course of the simulations we realized that there are appreciable sample to sample 
variations of the relaxation curves. Therefore it was necessary to average for each combination of Ti and Tf over 8-10 
independent runs. 



III. ONE-TIME QUANTITIES 

One-time quantities are observables which in equilibrium are constants, such as the total energy of the system, the 
pressure or, in a magnetic system, the magnetization. In the non-equilibrium situation these observables depend on 
the time which has elapsed since the quench and thus arc commonly called "one-time" quantities. 

It has been shown before, see e.g. Refs. [^jl5[l^, that the time dependence of such quantities is relatively weak. 



In Fig. |i| we show the time dependence of epot(i), the potential energy per particle, for various combinations of T 
and Tf. From the figure it becomes evident that this time dependence is rather weak. In Ref. |pl]] we have shown 
that it can be approximated well by a power-law with a small exponent (0.14), or alternatively by a logarithmic time 
dependence. Qualitatively this result agrees with the one of Monte Carlo simulations of a soft sphere system jlj] and 
a polymer model ]l5[ | in which also power-law dependences have been observed for one-time quantities. However, the 
exponents found in these simulations (0.7 and 1/3) are significantly larger than our, thus indicating that they are 
not universal quantities. The large exponent of reference |^2| might, however, be also due to the fact that in that 
simulation Ti was infinity and that the runs were rather short, i.e. that the dependence of Cpot on time was not the 
one valid for very long times. 

From Fig. pi we conclude that at long times the curves do not depend significantly on the value of the initial 
temperature (see the curves with Tf = 0.4). For short times, however, such a dependence can be seen, in that, e.g., 
the curve with T = 0.466 is almost constant and shows only at long times a time dependence which is similar to the 
one for higher values of Ti (see also Ref. [Q). This time dependence starts to be observable only for times that are 
comparable with the a-relaxation time of the system at Ti, which for Ti = 0.466 is on the order of 10** time units 
since this is the typical time scale that the system needs to find a significantly different configuration. This estimate 
is, however, only a lower bound, since the dynamics of the system is given by the temperature Tf < Ti, and is hence 
slower than the one at Ti. See also Ref. |l6| for a further discussion of this point. 

The figure also shows that there is a significant dependence of the curves on Tf in that the values of Cpot (t) decrease 
with decreasing Tf. In a harmonic solid at a temperature Tf one expects that the (equilibrium!) potential energy 
varies like 3/2kBTf and thus it is reasonable to subtract this "trivial contribution" from the curves. The result of this 
subtraction is shown in Fig. |l|b from which we see that this procedure does indeed make the curves collapse reasonably 
well. The most significant exceptions are the curves {Ti — 0.466, Tf — 0.4) and {Ti — 5.0, Tf = 0.1). These deviations 
can be understood by realizing that during the aging process the system is looking for configurations with low lying 
energies. If this search is started at a relatively low temperature, such as Ti = 0.466, the system will typically be 
able to find configurations that have a lower potential energy than in the case when Ti is large thus explaining why 
the curve for Ti = 0.466 is below the other ones. In addition to this, such a search will be more efficient at a higher 
temperature Tf than at a low one, since then the system has a better chance to overcome small local barriers. This 
explains why the curve for {Ti — 5.0, Tf = 0.1) is above the other ones with the same value of Ti. 

Since also the radial distribution function g{r) is a one-time quantity we expect that also its dependence on time is 
weak. That this is indeed the case is shown in the main part of Fig. y where we show gAA{T), i-e. the partial radial 
distribution function for the A particles (see Ref. ||l^ for its definition) for the times t — (i.e. before the quench), 
t = 10, 100, 1000, 10000, and 63100 time units. From the figure we see that immediately after the quench gAA{''') 
changes its shape very quickly in that, e.g., the first nearest neighbor peak which in the high temperature phase (curve 
labeled with < = 0) is not very high, becomes much higher and narrower. For times larger than 10, the shape becomes 
essentially independent of time, as expected for a one-time quantity and in agreement with the results of Ref. [p!2l. 
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That the form of gAA{f) o,t long times has a significant dependence on Tf is demonstrated in the inset of Fig. 
where we show the main peak for different values of Tf at long times {t = 63100). As expected we find that the height 
of this peak increases with decreasing temperature and that it becomes narrower. Thus we conclude that the typical 
configurations in which the system is stuck after long times does depend on the final temperature. Very recently Latz 
made the interesting proposition that the typical configurations in which the system is stuck at long times during 
the aging process share a common property in that all of them are very close to the so-called critical surface of 
mode-coupling theory (MCT) which divides the liquid like phase of the system from its glass like phase [^Sj. 
MCT predicts that this critical surface can be calculated by using as input the radial distribution function, i.e. a 
purely static quantity Furthermore it has been shown thatjfor simple liquids, the relevant part of this input is 

related to the area under the first peak in the structure factor , or alternatively, to the area under the first peak 
in g{r) weighted with Airr^ . In order to test the validity of the proposition by Latz we have therefore calculated the 
integral of 47rr^gAA{'') between zero and 1.406, the location of the first minimum in gK\(T). Note that this integral, 
which we call c, is the partial coordination number. The time and Tf dependence of c are shown in the main part 
of Fig. H (open symbols). (In order to expand the axis at low temperatures we show this data versus the logarithm 
of Tf.). The times t shown are spaced essentially equidistant on a logarithmic time axis and are 0, 10, 40, 60, 100, 
160, 250, 400, 630, 1000, 1580, 2510, 3980, 6310, 10000, 15850, 25120, 39810, and 63100. (In order to show the time 
dependence we plot the data for t = 63100 at T = Tf and with decreasing t at a temperature which is a factor of 1.003 
higher than the previous data point.) From this figure we recognize that with increasing time the value of c increases 
rapidly and then becomes constant to within the precision of our data. The value of this limiting constant, which we 
call Ccc, increases weakly with increasing Tf but to a first approximation it can be considered as independent of Tf, 
thus supporting the prediction of Latz. 

In order to see whether the Tf dependence of Coc is indeed weak it is useful to compare it with the temperature 
dependence of the area under the first peak ^aa in equilibrium. For this we have included in the figure also the value of 
c for temperatures 5.0 >T> 0.446 (filled symbols). (The gAAif) for the various temperatures are results from other 
simulations ^,|l^.) In this temperature range the equilibrium value of c shows an appreciable temperatures dependence 
(which is close to a logarithmic dependence) which is indeed much more pronounced than the Tf dependence of Coo of 
the non-equilibrium simulations. To a first approximation the value of Coo is given by the value of c of the equilibrium 
curve at T = Tc, where Tc is the so-called critical temperature of MCT ||l^. For the present system this critical 
temperature has been estimated to be around T = 0.435 [p|pC|]. Due to the problem to equilibrate the system at 
temperatures close to T^ it is currently not possible to determine the equilibrium value of c at T^ and therefore to 
compare this value with the Coo from the out-of-equilibrium simulations. However, since the equilibrium value of c is 
expected to show a smooth temperature dependence, it can be estimated quite reliably from the equilibrium data at 
a bit higher temperatures. In the inset of Fig. ^ we show an enlargement of the region which is marked by a box in 
the main figure and which encloses the temperature range around Tc. From this inset we see that the values of Coo 
for Tf = Tc — 0.435 and Tf = 0.4 are indeed very close to an extrapolation of the equilibrium curve to Tc. Hence 
we conclude that it is indeed possible to calculate from the equilibrium data with reasonable accuracy also certain 
quantities for the out-of-equilibrium situation. This point is discussed in more detail in Ref. |l^ . 

IV. TWO-TIMES QUANTITIES 
A. General features of two time correlation functions 

In equilibrium, time correlation functions between any two observables A and B, {A{t)B{0)) depend only on the 
time difference t, i.e. {A{t)B{0)) — {A{t +tiu) B (1^^)) , for all waiting times t^^. For the out of equilibrium situation this 
equality no longer holds and therefore such time correlation functions depend on two quantities, the time difference 
T and tyi,, the time passed since the generation of the non-equilibrium situation. Therefore such correlation functions 
are called two-times quantities and in this subsection we will demonstrate that such quantities are very well suited to 
investigate the aging dynamics of out of equilibrium systems. 

For liquids in equilibrium the dynamics is often studied by means of the intermediate scattering function F(k,t) 
which is defined by |l^ 

and by the its so-called self part, Fs{k,t) given by 
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Here rfc(T) are the positions of the particles at time t and k is the wave vector. For the non-equihbrium situation 
these functions are generaUzed to 

and 

j 

(Note that these last equations are trivially generalized to multi-component systems.) For the present system the 
dependence of F{k,t) and Fs{k,t) on time and temperature have been discussed extensively in Refs. [p|,p^. In Fig. ^ 
we show the time dependence of the corresponding quantities Ck{tw + t, tw) and Cc,k{tw + tw) for different waiting 
times t^,. The values of t^, are 0, i.e. immediately after the quench, and = 10, 100, 1000, 10000, and 63100 time 
units and the final temperature Tf = 0.4 (solid lines). In Fig. ||a we show Ck{tw + tw) for k — 7.23, the location of 
the first peak in the partial structure factor of the A-A correlation Q . We see that the correlation function decays 
quite quickly as a function of r, which in view of the fact that we are at a very low temperature might be surprising. 
However, one should realize that the decay of the correlation function to zero does not imply that the system has 
relaxed to equilibrium but only that the particle configuration at the end of the measurement is uncorrelated with 
the configuration at the start. From the figure we see that with increasing value of the waiting time the decay of 
the curves occurs at longer and longer times. This observation can be rationalized by realizing that the driving force 
which leads to the decorrelation of the state at time r from the state at time zero decreases with increasing waiting 
time since in the time between zero and t^j the system had already the possibility to relax. This dependence 
of the curves is, however, only observed at long times. For short times the curves fall onto a master curve and leave 
this master curve only at a time which increases with (and is roughly proportional to [Tl] ] ) . The time regime in 
which this master curve is observed is usually called "short time regime" whereas the time regime in which the curves 
show a significant 1^ dependence is called the "aging regime" . In the short time regime the particles rattle in the 
cages formed by their nearest neighbors and this type of motion is thus not sensitive on the value of (if is not 
too small). Only for longer times the particles are able to leave this cage and the typical time scale for this process 
depends strongly on f^. 

We also mention that the oscillation seen in the curves at r ~ 1.0 originates from the coupling of the system to 
the heat bath. Since at this time the velocities of all particles are swapped with the ones drawn from a Maxwell- 
Boltzmann distribution at temperature Tf, the motion of the particles is, on average, slowed down or even reversed. 
Hence for a brief period the relaxation is slower than expected and thus the correlation curve decays slower. This 
effect is, however, only seen at short times and thus can be disregarded for large values of r. 

It is also interesting to compare the relaxation behavior of the system in this non-equilibrium situation with the 
one in equilibrium. In Fig. ^ we have therefore also included the incoherent intermediate scattering function Fs{k, t) 
in equilibrium (dashed lines) which was obtained in a simulation of Gleim et al. ]T^ . The dashed line is a curve that 
corresponds to the equilibrium relaxation dynamics of the system at T = 0.446 (this was the lowest temperature at 
which the system could be equilibrated). For very short times, t < 1.0, this curve essentially coincides with the above 
discussed master curve in the short time regime. For times r larger than 1.0 deviations are observed, which are not 
due to the mentioned oscillations. A careful inspection of the curves shows that the approach of the master curve and 
the dashed curve to the (quasi) plateau at intermediate times is quite different from each other in that the former 
shows a very slow approach whereas the one of the latter is quite fast. 

Also at long times we find differences in the relaxation behavior. In equilibrium it has been demonstrated 
that the relaxation curves can be fitted very well by the so-called Kohlrausch- Williams- Watts (KWW) function, 
v4exp(— (t/ti-oi)^), where Tj^^i is the relaxation time and an exponent /? < 1 (sj. From the figure we see that at long 
times the shape of the equilibrium and non-equilibrium curves are quite different and a more detailed analysis shows 
that a KWW function does not give a good fit to the data. What seems to be the same, however, (or at least very 
similar) is the height of the plateau, i.e. the so-called Edwards- Anderson parameter ||2l[| or non-ergodicity parameter. 

In Fig. ^ we show the relaxation curves for the collective function, i.e. Cc,k{tw + t, t) for the same waiting times. 
The wave- vector is now k = 9.60, the location of the minimum in the partial structure factor of the A-A correlation 
Due to the collective nature of this correlation function the statistics is worse than the one for the self part but from 
the figure we see that also this observable can be studied with satisfactory accuracy. Since at this wave-vector the 
(equilibrium) value of the non-ergodicity parameter has a local minimum |^] it is expected that also the one for the 
non-equilibrium case is relatively small and from the figure we see that this is indeed the case. 
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B. Quantitative analysis of the relaxation at high temperatures 



In the following we will now analyze the relaxation behavior of the system in more detail. For this we start first 
with the dynamics at short times. Mean field calculations predict that for times at which the correlators are in the 
vicinity of the plateau, a time regime which in the equilibrium dynamics of glasses is called "/3-relaxation regime, two 
power-laws are observed [0,|2|. Denoting by qea the value of the Edwards- Anderson parameter, it is predicted that 
any time correlation function C(ttt, + r^t^) should behave like 

T if C > gEA (5) 

and 

C(t^, -I- r, t^) « gEA - Cb{T/Taf if C < ^ea ■ (6) 
Here is the typical relaxation time at long times and the exponents a and b are related by 

r2(i + 6) r2(i-a) 



T(l-h2&) r(l-2a)" 



(7) 



The quantity m in the last equation is the so-called fluctuation-dissipation-violation factor and will be discussed 
later in more detail. For the moment it is sufhcient to know that it is expected to be a system universal number 
equal or less than 1.0. We also mention that the two power-laws of Eqs. and (||) have been discussed for a long 
time and in great detail in the MCT of supercooled liquids and glasses pB. In order to see whether it is possible 
to see these power-laws in our data we have estimated the value of qea for Ck{tw + t) for k = 12.53 and show in 
Fig. ^ \Ck{tw + t) — (7ea(^)| in a double logarithmic plot for different waiting times. This value of k is relatively 
large so that the value of (/ea is relatively small, 0.47. A small value of (/ea is useful if one wants to check for the 
presence of the power-law given in Eq. (||) whereas for the check of Eq. (^) a large value of qea is better. (These 
conclusions are reached from the analysis of equilibrium dynamics but it can be expected that they hold also for 
the nonequilibrium case.) From the figure we recognize that it is indeed possible to see a power-law in the short time 
regime with an exponent around 0.42 ± 0.05 (bold solid line). In the discussion of Fig. ^ we said that in the late part 
of the /3-relaxation regime the equilibrium curves and the non-equilibrium curves show a very similar time dependence 



and that the one of equilibrium is given by a power-law. The exponent of this power-law is around 0.63 |14 2^]. If we 
assume the relation (^ to hold true one therefore obtains m « 0.57. Due to the relatively large error of the exponents 
a and h this number also has a significant error. 

We have also checked for the presence of the power-law for other values of fc, as well as for the collective correlation 
fimction Cc,fc(iu,-I-T, i^) and found that all of them show such a time dependence with an exponent which is compatible 
with 0.42. Unfortunately, however, the exact determination of the exponent is rather difficult, since a change of the 
(unknown) value of qea will lead to a change in the exponents also. In order to avoid this problem to some extend we 
use a trick which has proved to be useful in the context of the analysis of equilibrium data (see, e.g., p^ , p5|}_ . For the 
equilibrium case MCT predicts that in the ^-relaxation regime the so-called factorization property holds~]^,Q. This 
means that any time correlation function (/)(i) can be written as <j){t) = q-E,A{4') + h^G{t), where is a constant and 
the whole time dependence is given by the (/)-independent, i.e. system universal, function G{t). If we assume that a 
similar relation holds also in the non-equilibrium situation we have 

(l){tw+T,tw) =qEA{<P) + h^G{tiu+T,T) . (8) 

We see that if the exponents a and b as well as the quantity m are independent of the observable that then Eqs. (||) 
and (^ are indeed compatible with such an Ansatz. From Eq. (H) it follows immediately that if r' and r" are two 
arbitrary times in the /3-relaxation regime, the ratio 

n /, , X (t>(tw +T,t^) - (f){t^ + t", tt„) 

K4.[t^ + r, j - ^^^^ + _ ^(i^ + (9) 

is independent of if r is in the /3-relaxation regime also. In Fig. ^ we plot this ratio for different choices of 4>. 
These are Cc,k{tw + T^t^) for k ~ 6.52, 7.23, 9.6, and 12.53 as well as Ck{tw + T,tw) for the same wave-vectors 
and also k = 2.0 and k = 3.0. The value of is kept fixed at = 63100. From the figure we see that in the 
/3-relaxation regime the different curves collapse nicely onto a master curve, thus demonstrating the validity of the 
factorization property. That this result is by no means trivial can be recognized from the fact that for very short and 
very long times the different curves do not fall onto a master curve at all. From the existence of the master curve in 
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the /5- relaxation regime we thus come to the conclusion that, within the accuracy of our data, in this time regime the 
time dependence of the relaxation is governed by one single system universal function G{tw + t, twY- From the results 
shown in Fig. |^ we see that this function is compatible with power- law of the form given by Eq. (||) . 

As we will see later, it is expected that the quantity m depends on the value of Tf [0. Therefore one might conclude 
that also the exponent a and b depend on temperature. From making fits to the master curves in the /3-relaxation 
regime of our data, see Fig. it is hard to conclude whether or not such a dependence is indeed present since the 
choice of the non-ergodicity parameter affects the values of the exponents also. Therefore we have calculated the ratios 
R<f,{tw + T,tw) also for lower values of Tf, Tf = 0.3 and 0.1, and found that the curves for the different observables 
do indeed fall onto a master curve also. More important is the observation that within the accuracy of our data this 
master curve does not depend on Tf, thus giving evidence that the function G'(tu, + t, t,„) is only a weak function of 
Tf. ^ 

We also mention that mean-field calculations lead to the expectation that the quantities a, b and m in Eq. (|^) do 
depend on the final temperature. (This is in contrast to the equilibrium MCT where m = 1 and a and b are independent 



of temperature.) In Sec. V C we will show that m shows a significant dependence on the final temperature in that it 
becomes smaller with decreasing Tf. Thus that result might seem to be in contradiction with the fact that G(tw+T, t^) 
seems to depend only weakly on Tf. The resolution of this apparent problem is that if we assume b to be constant 
that a decrease of m can be compensated by an increase of a. It is, however, simple to see that if m is small a has to 
be close to 0.5 and that a variation of m can be compensated by a very small change in a, i.e. without leading to an 
appreciable change in G{t^ -|- r, 

Having analyzed the relaxation behavior of the system in the /3-relaxation regime we now investigate the one of the 
a-regime, i.e. the relaxation of the system at long times. When we discussed the time correlation functions in Fig. ^ 
we already mentioned that in this regime the relaxation differs from the one in the equilibrium situation in that the 
time correlation functions can not be described well by the KWW-law. In Fig. ^ we show the same correlation 
functions as in Fig. ^ja, but this time in a double log plot. From this presentation of the data we see that at long 
times the out-of-equilibrium curves show a power-law. The exponent is around 0.4 and is, within the accuracy of our 
data, independent of ^u, . A similar time dependence has also been found in simulations of spin-glasses p6|l , although 
in that case the exponent was significantly smaller, and is thus not unusual for the aging dynamics. 

The curves in Fig. ^ are for the wave-vector q = 7.23, the location of the main peak in the static structure factor. 
In Fig. ^ we show the same type of correlation function for different wave-vectors and = 1000. From the main 
figure we see that the dependence of Cfc on k is quite pronounced, in that e.g. the height of the plateau increases 
quickly with decreasing wave- vector. Such a dependence can be understood at least qualitatively by recalling that in 
equilibrium the wave- vector dependence of the Edwards- Anderson parameter is very similar to a Gaussian, and that 
at short times and long waiting times the time correlation functions do not depend on (see Fig. i.e. show a 
quasi-equilibrium behavior. From Fig. |^ we see that a rough estimate of qk is given by the value of Ck{tw +T,tw) 
at r = 20 and we find that this quantity does indeed show a Gaussian like dependence on k. 

In the inset of the figure we show the same correlation functions in a double logarithmic plot. From this graph we 
recognize that for all wave-vectors the long time dependence of the functions are compatible with a power-law and 
that the exponent depends on the wave-vector in that it decreases with decreasing k. 

Since at long times Ck shows a power-law dependence it is not possible to define a relaxation time. However, 
from the figure it becomes clear that the relaxation of the system is much faster on small length scales than on large 
ones. It is instructive to recall that for a diffusion process the relaxation time depend on the wave- vector like k~^. 
A comparison of the curves for k — 3.0 and k — 6.5 shows that the time needed to decay to 50% of the initial value 
differs by almost a factor of 100, thus much more than the factor of 4 expected for a diffusive process, or the factor 
found in the relaxation dynamics of supercooled liquids in equilibrium. This shows that during the aging process the 
relaxation at long times is indeed very different from the one in supercooled liquids. 

We now turn our attention to the time and waiting time dependence of the two-time correlation functions at long 
times. Within mean-field theory it is expected that for systems which show a discontinuous dependence of the non- 
ergodicity parameter on temperature, which is the case for the structural glass studied here, this dependence can be 
written (in the limit of long waiting times) as 

C{t^ + T,t^) = CsT(r) + Cag ( ^^^^^^^^ ) ■ (10) 

Here Cst(t) is the time dependence of C(it„ + T,t^) at short times which is supposed to be independent of (in 
agreement with our findings, see Fig. ^) and to decay quickly to zero. Cag is a function whose form depends on 
C{tw -\-T,tw) and h{t) is a monotonously increasing function of the argument. The interesting point of this equation 
is that the whole r and dependence of the aging regime enters only through the combination h{tw -\- T)/h{tw)- 
Apart from Eq. (Eo|), not much more is known about the t and dependence of C{t^ + r, t^). 
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In order to gain more insight into this dependence we investigate whether the function h{t) is independent of 
the correlation function C(t,i, + T,tyj). The presence of such a dependence can, e.g., be checked by plotting a time 
correlation function Ck'{tw + ''', ^w) versus a different time correlation function, e.g., Ck{tw +T,t^^), i.e. by making a 
parametric plot with time r as the running parameter. If h(t) is independent of the correlation function considered, 
it is easy to see that in such a parametric plot the curves corresponding to different waiting times will fall on top of 
each other. (This holds true only if the time dependence of Csrij) can be neglected, i.e. at long times.) 

In Fig. ^ we show such parametric plots for the waiting times tyj = 0, 10, 100, 1000, 10000, and 63100 time units. 
The independent variable (abscissa) is Ck{tw -'rT^tw) for k = 7.23 and the dependent variables (the ordinate) are the 
same correlation function for different values of k. Let us focus first on fc = 2.0. From the graph we recognize that 
at short times (corresponding to large values of and C^') the parametric plot is independent of t^, as it would 
be expected for the equilibrium case. However, for times at which the correlation functions have fallen below their 
Edwards- Anderson parameters a waiting time dependence is seen. Qualitatively the same behavior is found for the 
other values of k' . From this observation we hence conclude that if the Ansatz in Eq. (|lO|) is correct, then the function 
h{t) does depend on the observable, i.e. in this case on the wave-vector. We also mention that qualitatively the same 
results are obtained for the different values of T/. 

The time dependence of the function h{t) can be used to distinguish between different theoretical models to describe 
the aging dynamics. E.g. the droplet model of Fisher and Huse predicts h{t) to be of the form h{t) = log(t). 
In Ref. pS] it was argued that the present Lennard- Jones model showed this type of aging dynamics, a conclusion 
which was not confirmed by the present authors [p9| . The reason for this discrepancy remains still unresolved. In 
order to settle this issue we have attempted to determine the function h{t) from our data without making reference 
to any model, i.e. functional form of h{t). For this we assumed that at long times the correlation functions are given 
by the second term in Eq. Starting from the initial guess hit) = \og{t) we plotted the correlation function versus 
h{tw + T)/h{tw) for all values of i^,, which gave at long times a clustering of the curves. By iteratively bending h{t) 
in a continuous way by small amounts and subsequently making the scaling plot we attempted to generate a collapse 
of the curves. The outcome of this procedure is shown in Fig. ^ for the case Ck with k — 7.23. (For the sake of 
clarity only the curves for — 0, 10, 100, 1000, 10000, and 63100 time units are shown, although more waiting times 
were considered to determine the master curve. Also, in order to expand the abscissa we have subtracted 1.0 from 
h{tw + t) / h{tw) and plot its logarithm.). From this figure we see that the it is indeed possible to obtain a satisfactory 
scaling of the curves at long times. The function h(t) which was obtained by the procedure described above is shown 
in the inset. We see that to a first approximation it is close to a logarithm, but that significant deviations are present. 

Also included in the figure are the curves for different values of the wave- vector. We see that in these cases the 
curves for the different waiting times do not collapse nicely at long times. Thus this is further evidence that the 
function h{t) depends on the observable considered. We also mention that although we have determined the optimal 
shape of h also for the other wave-vectors, it is difficult to compare these different functions with each other. The 
problem is, that if h{t) is an optimal function that for example also hit)" (with arbitrary a) is an optimal function, 
i.e. there is no unique representation of h{t). 



C. Relaxation at low temperatures 



Essentially all the results discussed so far were obtained for a final temperature Tf = 0.4. This is only slightly 
(10%) less that the mode coupling critical temperature for the system, Tc = 0.435 ||]. A markedly different behavior 
in the relaxation is observed for temperatures much lower than Tc, as illustrated in figure |l^, which shows the results 
for Ck{tw -'r Zi'tw) after a quench to = 0.1. The values of the waiting times and of the wave- vector k are the same 
as in figure ^ {Tj = 0.4). At short times r, the dependences are qualitatively similar for the two final temperatures. 
The only difference is the increase in the plateau value when the temperature is lowered. This increase is easily 
understood from the fact that the amplitude of the vibrational motion about an equilibrium position decreases when 
T decreases. From harmonic theory a linear dependence of 1 — qsA on T at low temperatures can be expected, and 
such a dependence is indeed compatible with our results. 

For large values of and r, the behaviour for Tf = 0.1 is strikingly different from what was observed at Tf = 0.4. 
Instead of decaying rapidly to zero for t > t^, the correlators level off and appear to display an additional plateau at 
a value of C smaller than qEA- 

In order to get some insight into this surprising behavior, we have studied separately the correlation functions 



obtained for various samples prior to averaging. Looking at the data shown in fig. 10 for some of the samples, 
it appears that the decay to this second plateau is triggered by a large amplitude and rather sudden drop of the 
correlation function, that, depending on the sample, takes place typically 10'^ to 10'' time units after the quench (see 
inset of Fig. 11^). An analysis of the configurations shows that this sudden drop is related to a very collective dynamical 
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event, in which a substantial fraction of the particles (10%) move by a rather small amount, typically 0.1-Q.5crAA- 
Such motions can be understood from the fact that a deep quench leaves the system in a highly stressed configuration, 
so that the relaxation first proceeds through large "earthquake like" events, that release the local stress significantly. 
Only at longer times the aging dynamics crosses over to the very smooth process observed at higher values of Tf. 

We emphasize that these events do not seem to be related to the hopping like motion of the particles which is 
observed in the deeply supercooled state above the glass transition temperature fsot] . In that case only very few (2-3) 
particles are involved and during the jump these particles move on the order of one nearest neighbor distance. In 
contrast to this, the "catastrophic events" during the aging are very collective (50-100 particles) and the particles 
move only about 10-30 % of the typical nearest neighbor distance (as can be inferred from the self part of the van 
Hove correlation function). A typical event is shown in Fig. ^ where we show the particles just before the event 
(spheres) and their location after the event (tip of arrows). From this figure we recognize that in such an event the 
stress is released along a surface or even a line through the sample and not through the motion of the particles in a 
(bulk-like) three dimensional blob. 

Obviously, the catastrophic events that cause the decay of the correlations in this situation arc rather difficult 
to average over, and a very large number of samples would be required to obtain a reasonable statistical accuracy. 
Moreover, it is quite likely that these events correspond to a transient, so that the second plateau we observe is not 
really representative of the asymptotic behavior. Interestingly, however, rather similar shapes of the time correlation 
functions were observed by Bonn and coworkers | |3l[ | in their dynamical light scattering studies of aging in clay 
(laponite) suspensions. In that case the height of the second plateau was found to steadily increase with waiting time, 
confirming the transient nature of the effect. 

We also mention that the occurence of these events are related to the fact that the typical configuration at high 
temperature (T^ = 5.0) are quite different from the ones towards the system evolves to at Tf, thus giving rise to large 
stresses. Since in strong glass formers (e.g. Si02) the structure has a much weaker temperature dependance than in 
fragil eglassformers such as the present system, it can be expected that in the aging dynamics of the former no such 
"catastrophic events" should be observed. 



V. NONEQUILIBRIUM RESPONSE AND FLUCTUATION DISSIPATION RATIO 
A. Definition and measurement of the nonequilibrium response 

One of the crucial points that was derived from the solution of the dynamical equations in mean field models of 
spin glasses is the fact that, in a nonequilibrium system, the Fluctuation Dissipation Theorem (FDT) is violated in a 
nontrivial way. We will first briefly recall the main results from mean-field theories of spin glasses, and then discuss 
how the violation of FDT can be quantified in our system. 

Let us consider an observable A whose normalized autocorrelation function is denoted by C . If H denotes a field 
conjugate to A, the response oi A to H is defined as R{t,t') — j^jjr^- In an equilibrium system, this response 
function is invariant under time translation, i.e. R{t, t') = R{t — t'), and is related to the correlation function by the 
fluctuation dissipation theorem, R(t,t') = _I_ . Out of equilibrium, this relation does not hold any more, and 

a "Fluctuation Dissipation Ratio" (FDR) X{t,t') can be defined as 

RM^^xit,t'/-2Ml. (11, 

Much attention has been devoted to the asymptotic behavior of this fluctuation dissipation ratio. In particular, it 
has been shown |Q that in the asymptotic limit t,t' oo, the fluctuation dissipation ratio in mean field models of 
spin glasses, as well as in coarsening systems, becomes a non-singular function of C{t,t'), i.e. X{t,t') = x{C{t,t')). 
A direct consequence from this is that the more easily accessible integrated response 

M{t^+T,t„)^J R{t^+T,t)dt (12) 

can be written as a continuous function of C 

M{t,t')=M{C) [ x{c)dc . (13) 
Jc 
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From a practical point of view, the nontrivial consequence of these statements is that a parametric plot of M{tw + 
T,t^) versus C{tw + T,tw) (with t as the parameter) should, for long times, converge towards a master curve, 
independent of the waiting time. Such a behavior has been observed in a number of spin glass simulations p3|-05t. 

In order to test whether the same property holds in structural glasses, we have to devise a way of measuring the 
response function associated to Cfc(t, t'), which is the only correlation function that can be obtained with a reasonable 
accuracy in the aging regime. The procedure we use is the following. A Active "charge" e = ±1 is assigned randomly 
to each particle. An additional term of the form ejF(rj), where V^(r) — Vq cos(k-r) is a small (Vb < ksT) external 
potential, is then added to the Hamiltonian. It is then easy to check that, if one averages over several realizations of the 
random charge distribution, the time-correlation function of the observable Ak = exp(ik-rj(i)) is the incoherent 

scattering function C^. The procedure to generate the response function associated to is thus straightforward: For 
a given realization of the random charge distribution, the system is equilibrated at a high temperature (T!j = 5.0), 
and quenched at t = to the desired final temperature Tf. The evolution is followed with the field V{y) off until 
a waiting time t^, then the field is switched on and the response Ak{tw + T^tw) is monitored. The same procedure 
is repeated for several (7 to 10) realizations of the charge distribution, in order to get the response function. The 
quantity we obtain by this procedure is then an integrated response function M(tw + t, t„,) (Eq. since 

{Ak{t^ + T, ty,)) ^Vo R{tn, + T, t)dt (14) 

^VoM{t^+T,t^). (15) 

We have checked that this procedure indeed yields the correct response function by checking the fluctuation dissi- 
pation theorem in an equilibrium system. In that case, a slightly different procedure was used, in the sense that we 
monitored the decay of the response after switching off the field. The fluctuation dissipation theorem implies that 
this decay is directly proportional to the correlation function, a result that is illustrated in fig. |l^. From this figure 
we see that the quality of the response data is significantly poorer than the one for the correlation data. This is due 
to the fact that the latter can be averaged over time origins and wave-vector directions, hence improving drastically 
the statistical accuracy. 



B. Results for the response and the fluctuation dissipation ratio 

For the nonequilibrium case typical results for the integrated response function are shown in Fig. |l^, for Tf = 0.4, 
k = 7.25 and two values of the waiting time. As expected, the response, like the correlation, gets slower as the waiting 
time increases. As in the equilibrium case, the quality of the response data is poorer than that of the correlation 
data. In that case, the average over wave-vector directions in obtaining the correlation data probably explains this 
difference. Despite this noise it can be seen from the figure that at long times the curve for the response lies above 
the one for the correlation function, i.e. that the response is smaller than expected from the FDT. 

As explained above, mean field theories of spin glasses suggest that interesting information can be obtained from 
a parametric plot of M{tw + T,tw) versus C{tw + T,tw) at fixed tyj. Such a plot is shown in figures |lj (a-c) for 
three values of the final temperature, Tf = 0.4, Tf — 0.3 and Tf = 0.1. The wave-vector is fc = 7.25. For each 
temperature, results obtained at different waiting times are shown. In spite of the scatter in the data, the figures are 
clearly compatible with the existence of a master curve independent of the waiting time. Obviously, the quality of 
the data does not allow a quantitative analysis of this master curve. We can nevertheless argue that two different 
regimes can be distinguished. For high values of the correlation (Cfc close to unity), the parametric plot is essentially a 
straight line with slope —1. A slope of —1 would be observed in an equilibrium system obeying FDT. At lower values 
of the correlation, a clear deviation from this FDT slope is observed. Although other fits are certainly possible, we 
can describe the parametric plot as consisting of two straight lines, one with slope —1 and one with a negative slope 
— 1 < — m < 0. The reason for choosing such a description and the corresponding interpretation will be discussed in 
the next section. 

A similar analysis was performed for another value of the wave- vector {k = 3.0). Very similar results are obtained, 
and the slope of m of the non-FDT part in the parametric plot seems to be independent of the wave-vector. 

Before discussing the results, we mention that the results for this parametric plot are very sensitive to the linearity 
of the response. Hence we always checked carefully that our results were independent of the amplitude of the applied 
field. In practice, applying a stronger field does not affect the FDT part of the parametric plot, but tends to affect 
strongly the second part, in that it leads to a decrease in the parameter m that describes the non-FDT part of the 
curve. 
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C. Discussion 



The importance of the fluctuation dissipation ratio x{C) in glassy systems was first recognized in the context of the 
mean field theory of spin glasses [p2[ , when it appeared that this FDR was intimately related to the nature of ergodicity 



breaking in the system. The nature of this relationship was recently clarified by Franz et al. |36|, who showed that 
the asymptotic value of the FDR could be quite generally expressed in terms of the probability distribution P{q) of 
overlaps between replicas introduced by Parisi [p7| through 

x{C) = r P{q)dq. (16) 
Jo 

This relation between a purely static quantity, P{q), and a dynamical one, x{C), is extremely powerful, since 
it implies that a nontrivial information on the nature of phase space may be encoded in the nonequilibrium time 
dependent behavior. 



If we accept equation (16), a small number of possible scenarios arc documented in the literature |l38| . For systems 
with continuous step replica symmetry breaking, P{q) is a continuous function, and so is x{C). For systems with one 
step replica symmetry breaking like p-spins systems {p > 3), P{q) consists of two (5-functions at q = and q — qsA, 
so that x{C) = 1 for C > qsA and x{C) = m < 1 for C < qsA- Finally simple coarsening systems, like the Ising 
system below its ferromagnetic transition, have an essentially trivial P{q), P{q) — S{q — M^), where M — M{T) is 
the magnetization. Therefore, in that case the FDR is 1 if 1 > C > M^, and if > C |3|-0. 

Within this theoretical framework, we have to find which scenario is compatible with the results presented in the 
last section. From Figs. ^ it would appear that the most likely scenario is that of a system with one step replica 
symmetry breaking, for which x{C) is a stepwise constant function. In that case the parametric plot consists in two 
straight lines, one with slope —1 and one with slope — 1 < —m < (bold straight lines). 

Obviously the asymptotic nature of the results can be questioned. However, the fact that we obtained an essentially 
tw independent plot is already an indication that we are approaching the asymptotic limit. It could also be that 
preasymptotic effects cancel out when the parametric plot is used. Indeed, such a plot is even obtained for Tf = 0.1, 
where we have seen that strong "catastrophic" events influence the relaxation. Finally, we mention that in systems in 
which preasymptotic effects have been observed and studied (mostly domain growth models) |^,^ , they very clearly 
show up in the parametric plot as a t^^ dependence of the crossover region between the FDT part and the non-FDT 
region. The overall shape of the parametric plot is not affected much. Such a dependence could then explain the fact 
that in our data, the crossover between the FDT part and the non-FDT part takes place at a value of C smaller than 
the plateau value qEA in the correlation function, a feature which contradicts theoretical expectations. The tendency 
of the curves to "overshoot" in the crossover region, perceptible in Fig. is reminiscent of the observations made 
in reference [Q, and could also be a transient effect. In that paper it was argued that, e.g., in a coarsening system 
with some defects the de-pinning of a domain wall from a defect will give rise to a strong enhancement of the response 
just after the event. Since the "catastrophic" events correspond to a violent and, most likely, pre-asymptotic release 
of the stress, it is reasonable to assume that also in this case the response will be larger than expected. 

If we accept that our parametric plots correspond to the one step replica symmetry breaking case, we can extract 
from the data an estimate for the slope — m of the non-FDT part. The resulting slopes are given by: Tf = 0.4, 
m = 0.62 ± 0.05; Tf = 0.3, m = 0.45 ± 0.05; Tf = 0.1, m = 0.2 ± 0.1. As mentioned above, these results appear to 
be independent of the wave-vector. Within the accuracy of our data, they are compatible with a linear dependence 
of m on Tf, quite similar to that found by Parisi |l^ for a soft-sphere system. We also note that the value found for 
Tf — 0.4 is compatible with the result m = 0.57 that was extracted from a scaling analysis of the correlation functions 
alone. 

It is interesting to analyze our results for m using the "effective temperature" concept introduced in reference 
p2t . (For a different approach to the effective temperature see the papers of Nieuwenhuizen |4^ ) . In this approach, 
an effective "fluctuation dissipation temperature" T^ff is deflned as T^ff = T/m, where T is the actual external 
temperature. Crudely speaking, the "fast" degrees of freedom (those that correspond to the rapid decay of C to its 
plateau value) are at equilibrium with the thermostat at temperature T, while the slow degrees of freedom, which 
govern the aging behavior, are at equilibrium at a higher temperature. Then the concept of effective temperature can 
be thought as a rationalization of the older "flctive temperature" concept Q]. Within this interpretation, a linear 
dependence of m{Tf) on Tf corresponds to a constant effective temperature. It is therefore quite natural to expect 
flrst that m{Tf) should be equal to Tf/Tc, where Tc is the mode coupling critical temperature, since systems cooled 
below Tc fall out of equilibrium at Tc- This result was proposed by Parisi on the basis of his simulations in 
a soft sphere system. Our results do not corroborate such a simple interpretation. The "effective temperature" in 
our aging system is substantially larger than Tc, typically T^ff ~ 0.7. Our results here are similar those of Alvarez 
et al. p3| for the p-spin model in that these authors found for a temperature a bit below Tc a value of m which is 
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significantly smaller than 1. Indeed, the "ideal" result m{Tf) = Tf/Tc can only be expected to hold for systems that 
would be cooled infinitely slowly, so that they would remain in equilibrium down to Tc, and for which ergodicity 
restoring hopping processes can be completely neglected. For a system that is cooled with a finite cooling rate, it is 
not surprising to find an effective temperature which is above Tc- 



VI. CONCLUSIONS 



In this paper, we have presented a detailed numerical study of the out of equilibrium relaxation ("aging") of a 
simple glass forming system. The time scales we have investigated are those allowed by the current possibilities of 
Molecular Dynamics simulations, typically 10~^ — 10~^s. We will now try to briefly summarize the main conclusions 
that can be inferred from these numerical results. 

First, it appears that static quantities are, on these time scales, only very weakly dependent of time. By static 
(or "one time") quantities we mean, as usual, those quantities that can be obtained from the knowledge of a single 
configuration of the system. In particular, we have shown that the pair correlation functions (and hence all the 
derived quantities like the energy or the pressure) appear to equilibrate very quickly after the quench. These one- 
time quantities show some sensitivity to the quench history and to the final temperature and we find that to a first 
approximation they correspond to configurations which are very close to the critical surface of (equilibrium) mode- 
coupling theory. Hence we have evidence that for this model and within the time span of the simulation the system 
is not able to penetrate this surface and hence remains in the liquid like part of configuration space. 

In contrast to this weak dependence, the two-time correlation functions are very sensitive to the aging time. When 
the quench temperature is relatively high (i.e. still close to the mode-coupling critical temperature that was identified 
for our system), the waiting time dependence of these functions can be described in rather simple terms. Functions 
obtained for various can be rescaled on a master curve, and an analysis of this curve yields results that are 
compatible with the predictions of mean field ("mode-coupling") theories of the glass transition formulated for the 
"p-spin" spin glass models. 

At lower final temperatures, the behaviour of the two-time correlation functions is much more complex. The 
relaxation is largely dominated by "catastrophic events" that involve a small, but collective, displacement of a large 
number of particles. These infrequent and large events are very difficult to average over, so that the statistical quality 
of the data for the correlation functions is rather poor. Due to the limited time range of the simulations, we do not 
know whether these events are transients effects releasing some initial stresses due to the quench, or constitute a 
genuine characteristic of low temperature relaxation. 

The flexibility of MD simulations allows an independent measurement of the two-time response functions, which 
at equilibrium would be related to the correlation functions through the fluctuation dissipation theorem (FDT). 
A parametric representation of the response versus the correlation allows to clearly distinguish between a "quasi- 
equilibrium" region in which the FDT holds, and a second regime in which it is violated. An essential feature is that 
this parametric plot does not (or, at most, very weakly) depend on waiting time, so that the "FDT" and "non-FDT" 
regimes do not correspond to time windows but rather to correlation windows. This remarkable feature, predicted 
by mean field theories of spin glasses to hold in the asymptotic limit, is observed here for structural glasses at finite 
times. This is likely related to the formal similarities between the mode-coupling theories of structural glasses, which 
are known to describe quite well the equilibrium behaviour of our system, and these mean field models of spin glasses. 
As the response/correlation relationship observed for our model resembles the one observed in spin systems with one 
step replica symmetry breaking, one may speculate that the phase space structures in both systems are also similar. 
If this is the case, a quite appealing scenario for the glass transition can be devised, that reconciles some of the old 
ideas of the Adams-Gibbs approach with the more modern mode-coupling scenario, as discussed in detail in reference 

!§■ 

Unfortunately, we presently do not have at our disposal a quantitative theory of nonequilibrium behaviour in struc- 
tural glasses. The theoretical framework currently provided by mean-field theories of spin glasses can be considered 
as a schematic model, and hence does not allow comparisons beyond the qualitative level. At this level, we find that 
the general behaviour of our system is in quite good agreement with this schematic model. The agreement (and 
disagreements) are somewhat reminiscent of what is observed when comparing equilibrium data with predictions of 
schematic mode-coupling theories. 

Finally, it is natural to inquire about the relevance of a work on aging phenomena based on a method that is limited 
to time scales smaller than 10~^s to the aging effects observed on much longer time scales. The general picture we 
obtain provides, however, some support for a scenario |46| that may be of general relevance for much longer time 
scales. In this scenario, the critical temperature Tc of mode coupling theory is associated to the dynamical freezing 
temperature of mean field spin models. In a purely mean field (or, in other words, ideal mode-coupling) situation. 
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the system would be out of equilibrium at all temperatures below Tc- Of course, real systems can still be equilibrated 
below Tc, as "hopping processes" allow a relaxation that is not described by mode-coupling theory (at least in its 
ideal version). Hence in real experiments (or in our system, if we could allow for longer simulation times) these 
systems will show only interrupted aging, i.e. aging for short waiting times and no aging at long waiting times. 
However, we have shown that if the system was artificially driven into a nonequilibrium situation, its "short time 
aging" can be described reasonably well within the framework of mean- field/mode-coupling theories. We may here 
draw a parallel with the fact that, in equilibrium, below Tc, when the a (terminal) relaxation becomes dominated by 
hopping processes, the /3 (intermediate times) relaxation is still very well accounted for by ideal mode coupling theory 
It is then quite tempting to speculate that for any "fragile" system (i.e. a system in which hopping processes 
are not too strong, so that mode coupling theory accounts well for the dynamics close to Tc) quenched below its glass 
transition, the mean-field/mode-coupling description will be relevant for a large fraction of the aging process. This 
should be particularly relevant in systems where hopping processes are known to be weak, like colloidal suspensions. 
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FIG. 1. a) Time dependence of the potential energy for various combinations of Ti and T/. b) The same date as in a) but 
shifted by 3/2T/. 
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FIG. 2. Main figure: radial distribution function gAA(r) for times t = (before the quench) and t = 10, 100, 1000, 10000, 
and 63100. Ti = 5.0, T/ = 0.4. Inset: gAA{r) at long times for Ti = 5.0 and different values of T/. 
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FIG. 3. Main figure: Tf and time dependence of the area under the first peak in grAA(r) (open symbols). Closed symbols: 
Temperature dependence of this quantity in equilibrium. Inset: Enlargement of the region maked by a box in the main figure. 
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FIG. 4. Time dependence of the correlation functions Ck{tw + Tjtw), pannel (a), and Cc,k{tw + T,tu,), pannel (b), for the 
waiting times tw=0, 10, 100, 1000, 10000, and 63100. In pannel (a) we have also included an equilibrium curve at higher 
temperature (dashed line). See text for details. 




FIG. 5. Test of the presence of a power-law at short times in Ck{tw + T,tw)- The value of ^ea is 0.47. tuj=10, 100, 1000, 
10000, and 63100. The bold straight line has slope 0.42. 
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FIG. 6. Time dependence of the ratio R^{tw + t, tw) for different correlation functions (see text). The vertical arrows show 
the values of the times r' and t" used to calculate R^. 
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FIG. 7. a) Double logarithmic plot of Ck{tw + T,tw) for the aging dynamics (solid lines) and the equilibrium dynamics 
(dashed curve). The waiting times are: tw=0, 10, 100, 1000, 10000, and 63100. b) Wave-vector dependence of Ck{tw +T,tw) 
for the waiting time t^, = 1000. From top to bottom the values of k are 2.0, 3.0, 6.5, 7.23, 9.6, and 12.5. Inset: the same curves 
in a double logarithmic presentation. 




FIG. 8. Parametric plot between different correlation functions (see axis labels). The different curves for given wave- vector 
correspond to different values of the waiting time. T/ = 0.4. 
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h(t^+T)/h(tJ-1 

FIG. 9. Main figure: The two time correlation function as a function of the variable /i(r + tw)/h{tw) — 1 for different 
wave-vectors and waiting times (see text). The function h{t) was chosen to make the curves for k = 7.23 collapse at long times 
and its time dependence is shown in the inset. 




FIG. 10. Main figure: Time dependence of the correlation functions Ck{tw +t, tw), for the waiting times tw = 10, 100, 1000, 
10000, and 39810. T/ = 0.1, k = 7.23. Inset: The same correlation function for tw = 1000 for the individual runs. 
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FIG. 11. Snapshot of the configuration just before (t = 5070, spheres) and just after (t = 7650, tip of axrows) the large drop 
in the time correlation function. In this event Ck{tw + t, tw), k = 7.23, decayed from 0.79 to 0.52. T/ = 0.1. 




FIG. 12. Test of FDT at equilibrium. The dashed curve is the intermediate scattering function Fs (k, t) for k = 7.23 and the 
solid curve is the corrsponding response, i.e. 1 — kBTM{k,T). The response data was obtained by averaging over 14 different 
runs. 
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FIG. 13. Integrated response function M and correlation function for the out of equilibrium case. 
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FIG. 14. Parametric plots a) Tf = 0.4 b) T/ = 0.3 c)T/ = 0.1 
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